module initSubgridMod

  !-----------------------------------------------------------------------
  ! !DESCRIPTION:
  ! Lower-level routines for initializing the subgrid structure. This module is shared
  ! between both the production code (via initGridCellsMod) and unit testing code.
  !
  ! !USES:
  use shr_kind_mod   , only : r8 => shr_kind_r8
  use shr_log_mod    , only : errMsg => shr_log_errMsg
  use spmdMod        , only : masterproc
  use abortutils     , only : endrun
  use clm_varctl     , only : iulog, use_fates
  use decompMod      , only : bounds_type, subgrid_level_landunit, subgrid_level_column, subgrid_level_patch
  use GridcellType   , only : grc                
  use LandunitType   , only : lun                
  use ColumnType     , only : col                
  use PatchType      , only : patch                
  use column_varcon  , only : is_hydrologically_active
  !
  ! !PUBLIC TYPES:
  implicit none
  private
  save
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  public :: clm_ptrs_compdown ! fill in data pointing down
  public :: clm_ptrs_check    ! checks and writes out a summary of subgrid data
  public :: add_landunit      ! add an entry in the landunit-level arrays
  public :: add_column        ! add an entry in the column-level arrays
  public :: add_patch         ! add an entry in the patch-level arrays
  !

  character(len=*), parameter, private :: sourcefile = &
       __FILE__
  !-----------------------------------------------------------------------

contains
  
  !------------------------------------------------------------------------------
  subroutine clm_ptrs_compdown(bounds)
    !
    ! !DESCRIPTION:
    ! Assumes the part of the subgrid pointing up has been set.  Fills 
    ! in the data pointing down.  Up is p_c, p_l, p_g, c_l, c_g, and l_g.
    !
    ! This algorithm assumes all indices besides grid cell are monotonically
    ! increasing.  (Note that grid cell index is NOT monotonically increasing,
    ! hence we cannot set initial & final indices at the grid cell level - 
    ! grc%luni, grc%lunf, etc.)
    !
    ! Algorithm works as follows.  The p, c, and l loops march through
    ! the full arrays (nump, numc, and numl) checking the "up" indexes.
    ! As soon as the "up" index of the current (p,c,l) cell changes relative 
    ! to the previous (p,c,l) cell, the *i array will be set to point down 
    ! to that cell.  The *f array follows the same logic, so it's always the
    ! last "up" index from the previous cell when an "up" index changes.
    !
    ! For example, a case where p_c(1:4) = 1 and p_c(5:12) = 2.  This 
    ! subroutine will set c_pi(1) = 1, c_pf(1) = 4, c_pi(2) = 5, c_pf(2) = 12.
    !
    ! !USES
    use clm_varcon, only : ispval
    !
    ! !ARGUMENTS
    implicit none
    type(bounds_type), intent(in) :: bounds  ! bounds
    !
    ! !LOCAL VARIABLES:
    integer :: l,c,p               ! loop counters
    integer :: curg,curl,curc,curp ! tracks g,l,c,p indexes in arrays
    integer :: ltype               ! landunit type
    !------------------------------------------------------------------------------

    !--- Set the current c,l (curc, curl) to zero for initialization,
    !---   these indices track the current "up" index.
    !--- Take advantage of locality of l/c/p cells
    !--- Loop p through full local begp:endp length
    !--- Separately check the p_c, p_l, and p_g indexes for a change in
    !---   the "up" index.
    !--- If there is a change, verify that the current c,l,g is within the 
    !---   valid range, and set c_pi, l_pi, or g_pi to that current c,l,g
    !--- Constantly update the c_pf, l_pf, and g_pf array.  When the
    !---   g, l, c index changes, the *_pf array will be set correctly
    !--- Do the same for cols setting c_li, c_gi, c_lf, c_gf and
    !---   lunits setting l_gi, l_gf.

    curc = 0
    curl = 0
    do p = bounds%begp,bounds%endp
       if (patch%column(p) /= curc) then
          curc = patch%column(p)
          if (curc < bounds%begc .or. curc > bounds%endc) then
             write(iulog,*) 'clm_ptrs_compdown ERROR: pcolumn ',p,curc,bounds%begc,bounds%endc
             call endrun(subgrid_index=p, subgrid_level=subgrid_level_patch, msg=errMsg(sourcefile, __LINE__))
          endif
          col%patchi(curc) = p
       endif
       col%patchf(curc) = p
       col%npatches(curc) = col%patchf(curc) - col%patchi(curc) + 1
       if (patch%landunit(p) /= curl) then
          curl = patch%landunit(p)
          if (curl < bounds%begl .or. curl > bounds%endl) then
             write(iulog,*) 'clm_ptrs_compdown ERROR: plandunit ',p,curl,bounds%begl,bounds%endl
             call endrun(subgrid_index=p, subgrid_level=subgrid_level_patch, msg=errMsg(sourcefile, __LINE__))
          endif
          lun%patchi(curl) = p
       endif
       lun%patchf(curl) = p
       lun%npatches(curl) = lun%patchf(curl) - lun%patchi(curl) + 1
    enddo

    curl = 0
    do c = bounds%begc,bounds%endc
       if (col%landunit(c) /= curl) then
          curl = col%landunit(c)
          if (curl < bounds%begl .or. curl > bounds%endl) then
             write(iulog,*) 'clm_ptrs_compdown ERROR: clandunit ',c,curl,bounds%begl,bounds%endl
             call endrun(subgrid_index=c, subgrid_level=subgrid_level_column, msg=errMsg(sourcefile, __LINE__))
          endif
          lun%coli(curl) = c
       endif
       lun%colf(curl) = c
       lun%ncolumns(curl) = lun%colf(curl) - lun%coli(curl) + 1
    enddo

    ! Determine landunit_indices: indices into landunit-level arrays for each grid cell.
    ! Note that landunits not present in a given grid cell are set to ispval.
    grc%landunit_indices(:,bounds%begg:bounds%endg) = ispval
    do l = bounds%begl,bounds%endl
       ltype = lun%itype(l)
       curg = lun%gridcell(l)
       if (curg < bounds%begg .or. curg > bounds%endg) then
          write(iulog,*) 'clm_ptrs_compdown ERROR: landunit_indices ', l,curg,bounds%begg,bounds%endg
          call endrun(subgrid_index=l, subgrid_level=subgrid_level_landunit, msg=errMsg(sourcefile, __LINE__))
       end if

       if (grc%landunit_indices(ltype, curg) == ispval) then
          grc%landunit_indices(ltype, curg) = l
       else
          write(iulog,*) 'clm_ptrs_compdown ERROR: This landunit type has already been set for this gridcell'
          write(iulog,*) 'l, ltype, curg = ', l, ltype, curg
          call endrun(subgrid_index=l, subgrid_level=subgrid_level_landunit, msg=errMsg(sourcefile, __LINE__))
       end if
    end do

  end subroutine clm_ptrs_compdown

  !------------------------------------------------------------------------------
  subroutine clm_ptrs_check(bounds)
    !
    ! !DESCRIPTION:
    ! Checks subgrid data
    !
    ! !USES
    use clm_varcon, only : ispval
    use landunit_varcon, only : max_lunit
    !
    ! !ARGUMENTS
    implicit none
    type(bounds_type), intent(in) :: bounds
    !
    ! !LOCAL VARIABLES:
    integer :: g,l,c,p       ! loop counters
    integer :: l_prev        ! l value of previous point
    integer :: ltype         ! landunit type
    logical :: error         ! error flag
    !------------------------------------------------------------------------------

    associate( &
         begg => bounds%begg, &
         endg => bounds%endg, &
         begl => bounds%begl, &
         endl => bounds%endl, &
         begc => bounds%begc, &
         endc => bounds%endc, &
         begp => bounds%begp, &
         endp => bounds%endp  &
         )
    
    !--- check index ranges ---
    error = .false.
    do g = begg, endg
       do ltype = 1, max_lunit
          l = grc%landunit_indices(ltype, g)
          if (l /= ispval) then
             if (l < begl .or. l > endl) error = .true.
          end if
       end do
    end do
    if (error) then
       call endrun( &
            msg = 'clm_ptrs_check: g index ranges - ERROR', &
            additional_msg = errMsg(sourcefile, __LINE__))
    end if

    error = .false.
    if (minval(lun%gridcell(begl:endl)) < begg .or. maxval(lun%gridcell(begl:endl)) > endg) error=.true.
    if (minval(lun%coli(begl:endl)) < begc .or. maxval(lun%coli(begl:endl)) > endc) error=.true.
    if (minval(lun%colf(begl:endl)) < begc .or. maxval(lun%colf(begl:endl)) > endc) error=.true.
    if (minval(lun%patchi(begl:endl)) < begp .or. maxval(lun%patchi(begl:endl)) > endp) error=.true.
    if (minval(lun%patchf(begl:endl)) < begp .or. maxval(lun%patchf(begl:endl)) > endp) error=.true.
    if (error) then
       call endrun( &
            msg = 'clm_ptrs_check: l index ranges - ERROR', &
            additional_msg = errMsg(sourcefile, __LINE__))
    endif

    error = .false.
    if (minval(col%gridcell(begc:endc)) < begg .or. maxval(col%gridcell(begc:endc)) > endg) error=.true.
    if (minval(col%landunit(begc:endc)) < begl .or. maxval(col%landunit(begc:endc)) > endl) error=.true.
    if (minval(col%patchi(begc:endc)) < begp .or. maxval(col%patchi(begc:endc)) > endp) error=.true.
    if (minval(col%patchf(begc:endc)) < begp .or. maxval(col%patchf(begc:endc)) > endp) error=.true.
    if (error) then
       call endrun( &
            msg = 'clm_ptrs_check: c index ranges - ERROR', &
            additional_msg = errMsg(sourcefile, __LINE__))
    endif

    error = .false.
    if (minval(patch%gridcell(begp:endp)) < begg .or. maxval(patch%gridcell(begp:endp)) > endg) error=.true.
    if (minval(patch%landunit(begp:endp)) < begl .or. maxval(patch%landunit(begp:endp)) > endl) error=.true.
    if (minval(patch%column(begp:endp)) < begc .or. maxval(patch%column(begp:endp)) > endc) error=.true.
    if (error) then
       call endrun( &
            msg = 'clm_ptrs_check: p index ranges - ERROR', &
            additional_msg = errMsg(sourcefile, __LINE__))
    endif

    !--- check that indices in arrays are monotonically increasing ---
    error = .false.
    do l=begl+1,endl
      if ((lun%itype(l) == lun%itype(l-1)) .and. &
           lun%gridcell(l) < lun%gridcell(l-1)) then
         ! grid cell indices should be monotonically increasing for a given landunit type
         error = .true.
      end if
      if (lun%coli(l) < lun%coli(l-1)) error = .true.
      if (lun%colf(l) < lun%colf(l-1)) error = .true.
      if (lun%patchi(l) < lun%patchi(l-1)) error = .true.
      if (lun%patchf(l) < lun%patchf(l-1)) error = .true.
      if (error) then
         call endrun(subgrid_index=l, subgrid_level=subgrid_level_landunit, &
              msg = 'clm_ptrs_check: l mono increasing - ERROR', &
              additional_msg = errMsg(sourcefile, __LINE__))
      endif
    enddo

    error = .false.
    do c=begc+1,endc
      l = col%landunit(c)
      l_prev = col%landunit(c-1)
      if ((lun%itype(l) == lun%itype(l_prev)) .and. &
           col%gridcell(c) < col%gridcell(c-1)) then
         ! grid cell indices should be monotonically increasing for a given landunit type
         error = .true.
      end if
      if (col%landunit(c) < col%landunit(c-1)) error = .true.
      if (col%patchi(c) < col%patchi(c-1)) error = .true.
      if (col%patchf(c) < col%patchf(c-1)) error = .true.
      if (error) then
         call endrun(subgrid_index=c, subgrid_level=subgrid_level_column, &
              msg = 'clm_ptrs_check: c mono increasing - ERROR', &
              additional_msg = errMsg(sourcefile, __LINE__))
      endif
    enddo

    error = .false.
    do p=begp+1,endp
      l = patch%landunit(p)
      l_prev = patch%landunit(p-1)
      if ((lun%itype(l) == lun%itype(l_prev)) .and. &
           patch%gridcell(p) < patch%gridcell(p-1)) then
         ! grid cell indices should be monotonically increasing for a given landunit type
         error = .true.
      end if
      if (patch%landunit(p) < patch%landunit(p-1)) error = .true.
      if (patch%column  (p) < patch%column  (p-1)) error = .true.
      if (error) then
         call endrun(subgrid_index=p, subgrid_level=subgrid_level_patch, &
              msg = 'clm_ptrs_check: p mono increasing - ERROR', &
              additional_msg = errMsg(sourcefile, __LINE__))
      endif
    enddo

    !--- check that the tree is internally consistent ---
    error = .false.
    do g = begg, endg
       do ltype = 1, max_lunit
          l = grc%landunit_indices(ltype, g)

          ! skip l == ispval, which implies that this landunit type doesn't exist on this grid cell
          if (l /= ispval) then
             if (lun%itype(l) /= ltype) error = .true.
             if (lun%gridcell(l) /= g) error = .true.
             if (error) then
                call endrun(subgrid_index=l, subgrid_level=subgrid_level_landunit, &
                     msg = 'clm_ptrs_check: tree consistent - ERROR', &
                     additional_msg = errMsg(sourcefile, __LINE__))
             endif
             do c = lun%coli(l),lun%colf(l)
                if (col%gridcell(c) /= g) error = .true.
                if (col%landunit(c) /= l) error = .true.
                if (error) then
                   call endrun(subgrid_index=c, subgrid_level=subgrid_level_column, &
                        msg = 'clm_ptrs_check: tree consistent - ERROR', &
                        additional_msg = errMsg(sourcefile, __LINE__))
                endif
                do p = col%patchi(c),col%patchf(c)
                   if (patch%gridcell(p) /= g) error = .true.
                   if (patch%landunit(p) /= l) error = .true.
                   if (patch%column(p)   /= c) error = .true.
                   if (error) then
                      call endrun(subgrid_index=p, subgrid_level=subgrid_level_patch, &
                           msg = 'clm_ptrs_check: tree consistent - ERROR', &
                           additional_msg = errMsg(sourcefile, __LINE__))
                   endif
                enddo  ! p
             enddo  ! c
          end if  ! l /= ispval
       enddo  ! ltype
    enddo  ! g

    end associate

  end subroutine clm_ptrs_check

  !-----------------------------------------------------------------------
  subroutine add_landunit(li, gi, ltype, wtgcell)
    !
    ! !DESCRIPTION:
    ! Add an entry in the landunit-level arrays. li gives the index of the last landunit
    ! added; the new landunit is added at li+1, and the li argument is incremented
    ! accordingly.
    !
    ! !USES:
    use landunit_varcon , only : istice, istdlak, isturb_MIN, isturb_MAX, landunit_is_special
    !
    ! !ARGUMENTS:
    integer  , intent(inout) :: li      ! input value is index of last landunit added; output value is index of this newly-added landunit
    integer  , intent(in)    :: gi      ! grid cell index on which this landunit should be placed
    integer  , intent(in)    :: ltype   ! landunit type
    real(r8) , intent(in)    :: wtgcell ! weight of the landunit relative to the grid cell
    !
    ! !LOCAL VARIABLES:
    
    character(len=*), parameter :: subname = 'add_landunit'
    !-----------------------------------------------------------------------
    
    li = li + 1

    lun%gridcell(li) = gi
    lun%wtgcell(li) = wtgcell
    lun%itype(li) = ltype

    lun%ifspecial(li) = landunit_is_special(ltype)

    if (ltype == istice) then
       lun%glcpoi(li) = .true.
    else
       lun%glcpoi(li) = .false.
    end if

    if (ltype == istdlak) then
       lun%lakpoi(li) = .true.
    else
       lun%lakpoi(li) = .false.
    end if

    if (ltype >= isturb_MIN .and. ltype <= isturb_MAX) then
       lun%urbpoi(li) = .true.
    else
       lun%urbpoi(li) = .false.
    end if

  end subroutine add_landunit

  !-----------------------------------------------------------------------
  subroutine add_column(ci, li, ctype, wtlunit, type_is_dynamic)
    !
    ! !DESCRIPTION:
    ! Add an entry in the column-level arrays. ci gives the index of the last column
    ! added; the new column is added at ci+1, and the ci argument is incremented
    ! accordingly.
    !
    ! !ARGUMENTS:
    integer  , intent(inout) :: ci      ! input value is index of last column added; output value is index of this newly-added column
    integer  , intent(in)    :: li      ! landunit index on which this column should be placed (assumes this landunit has already been created)
    integer  , intent(in)    :: ctype   ! column type
    real(r8) , intent(in)    :: wtlunit ! weight of the column relative to the landunit

    ! whether this column's type can change at runtime; if not provided, assumed to be false
    logical  , intent(in), optional :: type_is_dynamic
    !
    ! !LOCAL VARIABLES:
    logical :: l_type_is_dynamic  ! local version of type_is_dynamic

    character(len=*), parameter :: subname = 'add_column'
    !-----------------------------------------------------------------------

    l_type_is_dynamic = .false.
    if (present(type_is_dynamic)) then
       l_type_is_dynamic = type_is_dynamic
    end if

    ci = ci + 1

    col%landunit(ci) = li
    col%gridcell(ci) = lun%gridcell(li)
    col%wtlunit(ci) = wtlunit
    col%itype(ci) = ctype
    col%lun_itype(ci) = lun%itype(li)
    col%type_is_dynamic(ci) = l_type_is_dynamic
    col%hydrologically_active(ci) = is_hydrologically_active( &
         col_itype = ctype, &
         lun_itype = lun%itype(li))
    col%urbpoi(ci) = lun%urbpoi(li)

  end subroutine add_column

  !-----------------------------------------------------------------------
  subroutine add_patch(pi, ci, ptype, wtcol)
    !
    ! !DESCRIPTION:
    ! Add an entry in the patch-level arrays. pi gives the index of the last patch added; the
    ! new patch is added at pi+1, and the pi argument is incremented accordingly.
    !
    ! !USES:
    use clm_varcon      , only : ispval
    use landunit_varcon , only : istsoil, istcrop
    use clm_varpar      , only : natpft_lb
    !
    ! !ARGUMENTS:
    integer  , intent(inout) :: pi    ! input value is index of last patch added; output value is index of this newly-added patch
    integer  , intent(in)    :: ci    ! column index on which this patch should be placed (assumes this column has already been created)
    integer  , intent(in)    :: ptype ! patch type
    real(r8) , intent(in)    :: wtcol ! weight of the patch relative to the column
    !
    ! !LOCAL VARIABLES:
    integer :: li        ! landunit index
    integer :: lb_offset ! offset between natpft_lb and 1
    
    character(len=*), parameter :: subname = 'add_patch'
    !-----------------------------------------------------------------------
    
    pi = pi + 1

    patch%column(pi) = ci
    li = col%landunit(ci)
    patch%landunit(pi) = li
    patch%gridcell(pi) = col%gridcell(ci)

    patch%wtcol(pi) = wtcol

    ! TODO (MV, 10-17-14): The following must be commented out because
    ! currently patch%itype is used in CanopyTemperatureMod to calculate 
    ! z0m(p) and displa(p) - and is still called even when fates is on

    !if (.not. use_fates) then
    patch%itype(pi) = ptype
    !end if

    if (lun%itype(li) == istsoil .or. lun%itype(li) == istcrop) then
       lb_offset = 1 - natpft_lb
       patch%mxy(pi) = ptype + lb_offset
    else
       patch%mxy(pi) = ispval
    end if
    

  end subroutine add_patch


end module initSubgridMod
